home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / procssng / ccs / ccs-11tl.lha / lbl / hips / sources / 3d-tools / 3dslice.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-06-15  |  12.2 KB  |  454 lines
  1. /*
  2. % 3DSLICE.C - slicing 3D image to 2D image for displaying and looking the
  3. %    inside of the image. Default is to slice horizontally (x-z surface,
  4. %    parallel X axle). -dv specify slicing vertically (y-z surface, Z axle).
  5. %    -dd doing slicing diagnally and -p XYZ gives choice for parallel to
  6. %    which axle, and -a # gives angle which range is -90 -- +90 degree.
  7. %    But you can specify any angle, the program can convert it into this
  8. %    range.
  9. %
  10. %    Copyright (c)    1990    Jin, Guojun
  11. */
  12. char    usage[]="options    \n\
  13. -a#        slice angle    \n\
  14. -d[h] [v] [d]    slice direction    \n\
  15. -p[x] [y] [z]    parallel to axle\n\
  16. -c#        begin column    \n\
  17. -r#        begin row    \n\
  18. -f#        begin frame    \n\
  19. -n#        number of frames will be processed\n\
  20. input [> output]\n\
  21. % Note:    \n\
  22. %    The result of slicing is always from bottom to top.    \n\
  23. %    To expain this, we use following example:        \n\
  24. %    Parallel Y axle slicing. Coordinate is (x, y, z);    \n\
  25. %    (0, 0, 0 - 0, 1, 0) corner by using +angle -f options.    \n\
  26. %    (1, 0, 0 - 1, 1, 0) corner by using -angle -c options.    \n\
  27. %    (1, 0, 1 - 1, 1, 1) using +angle & -c options.        \n\
  28. %    (0, 0, 1 - 0, 1, 1) using -angle & -f options.        \n\
  29. %    For 1 & 3, the view-line is (1, 0, 1 - 1, 1, 1).    \n\
  30. %    For 2 & 4, the view_line is (1, 0, 0 - 1, 1, 0).    \n\
  31. %    The other surfaces will be similar.    \n";
  32. /*
  33. * AUTHOR:    Guojun Jin - LBL    11/15/90
  34. */
  35.  
  36. #include "header.def"
  37. #include "imagedef.h"
  38. #include <math.h>
  39.  
  40. U_IMAGE    uimg;
  41.  
  42. #define    inbuf    uimg.src
  43. #define    obuf    uimg.dest
  44. #define    cols    uimg.width
  45. #define    rows    uimg.height
  46.  
  47. double    ssin, sinl, cosl, tanl;
  48.  
  49. #ifndef    GValue
  50. #define    GValue()    arget(argc, argv, &i, &c)
  51. #endif
  52.  
  53.  
  54. main(argc, argv)
  55. int    argc;
  56. char**    argv;
  57. {
  58. byte    *inbp, *obp, *tmpp;
  59. bool    bc=0, bf=0;    /* default = row    */
  60. int    *ibuf, *ibp,    /* for diagnal interpolation    */
  61.     slice_dir, l_r=1,/* default for parallel to Z axle    */
  62.     xyz='z', width,
  63.     bgn_row=1,    /* all working number is start from 1    */
  64.     bgn_cln=0,    /* and keeps in all diagnal routines    */
  65.     bgn_frm=0,    /* and changed to from 0 at right angle routine    */
  66.     f, r, c, df,
  67.     rowsNew, clnsNew, fsizeNew;
  68. MType    i, fsize, frames=1;
  69. float    angle=45, rfmod, cfmod;
  70.  
  71. format_init(&uimg, IMAGE_INIT_TYPE, HIPS, HIPS, *argv, "S20-1");
  72.  
  73. for (i=1; i<argc; i++)
  74.     if (*argv[i] == '-' || *argv[i] == '+') {
  75.     c = 1;
  76.     switch (argv[i][c++]){
  77.     case 'a':
  78.         if (avset(argc,argv,&i,&c,0))
  79.             angle = atof(argv[i]+c);    break;
  80. #ifdef    _DEBUG_
  81.     case 'D':    debug++;    break;
  82. #endif
  83.     case 'c':
  84.         bc = bgn_cln = GValue();
  85.         bgn_frm = 0;    break;
  86.     case 'f':
  87.         bf = bgn_frm = GValue();
  88.         bc = 0;    break;
  89.     case 'd':
  90.         slice_dir = argv[i][c];    break;
  91.     case 'n':
  92.         frames = GValue();    break;
  93.     case 'p':
  94.         slice_dir = xyz = argv[i][c];    break;
  95.     case 'r':
  96.         bgn_row = GValue();
  97.         bf=0;    break;
  98.     case 'o':
  99.         if (out_fp=freopen(argv[i]+2, "wb", stdout))    break;
  100.     default:
  101. errout:        usage_n_options(usage, i, argv[i]);
  102.     }
  103.     }
  104.     else if ((in_fp = freopen(argv[i], "rb", stdin)) == NULL)
  105.         syserr("input file %s not found", argv[i]);
  106.  
  107. io_test(fileno(in_fp), goto    errout);
  108.  
  109. angle -= (int)(angle/180) * 180;
  110. if (fabs(angle) > 90)
  111.     angle -= 180;
  112. if (!angle)     slice_dir = 'h';
  113.  
  114. (*uimg.header_handle)(HEADER_READ, &uimg, 0, 0);
  115. if (uimg.in_form != IFMT_BYTE)
  116.     syserr("Non BYTE format isn't available Now! Use 3drotate please");
  117.  
  118. uimg.pxl_out = uimg.pxl_in;
  119. uimg.o_form = uimg.in_form;
  120. fsize = cols * rows;
  121.  
  122. inbuf = nzalloc((MType)frms, fsize, "inbuf");
  123. /* for right angle only, diagnal routines will alloc obuf by themself    */
  124. if (slice_dir=='h' || slice_dir=='v' || !slice_dir)
  125.     obuf = nzalloc(frames, fsize, "obuf");
  126.  
  127. /* parallel to Z axle only needs 1 integer bufffer,
  128.     and parallel to X or Y needs 2 buffers    */
  129. ibuf = (int*)zalloc(fsize, (MType)(sizeof(*ibuf) << (xyz=='x' || xyz=='y')), "ibuf");
  130.  
  131. (*uimg.std_swif)(FI_LOAD_FILE, &uimg, uimg.load_all=frms, No);
  132.  
  133. switch(slice_dir)
  134. {
  135. case 'v':
  136.     if (--bgn_cln < 0)    bgn_cln = cols/3;
  137.     if (frames > cols-bgn_cln)    frames = cols - bgn_cln;
  138.     i = cols;    cols = frms;
  139.     frms=frames;    fsizeNew = rows*cols;
  140.     (*uimg.header_handle)(HEADER_WRITE, &uimg, argc, argv, True);
  141.  
  142.     message("vertical at %d column\n", bgn_cln);
  143.     for (c=0; c<frames; c++)
  144.     {
  145.         tmpp = (byte*)inbuf + c + bgn_cln;
  146.         obp = (byte*)obuf;
  147.         for (r=0; r < rows; r++, tmpp+=i)
  148.             for (f=0, inbp=tmpp; f<cols; f++, inbp += fsize)
  149.             *obp++ = *inbp;
  150.         if(fwrite(obuf, uimg.pxl_out, fsizeNew, out_fp) != fsizeNew)
  151.             syserr("error during write.");
  152.     }
  153.     break;
  154. case 'h':
  155. case 0:
  156.     if (--bgn_row < 0)    bgn_row = rows/3;
  157.     if (frames > rows-bgn_row)    frames = rows - bgn_row;
  158.     rows = frms;
  159.     frms=frames;    fsizeNew = rows*cols;
  160.     (*uimg.header_handle)(HEADER_WRITE, &uimg, argc, argv, True);
  161.  
  162.     message("horizontal at %d rows\n", bgn_row);
  163.     for (r=0; r<frames; r++)
  164.     {    inbp = (byte*)inbuf + (r+bgn_row)*cols + fsize*(rows-1);/* rows = old_img.frames */
  165.         obp = (byte*)obuf;
  166.         for (f=0; f < rows; f++, inbp -= fsize + cols)
  167.             for (c=0; c < cols; c++)
  168.             *obp++ = *inbp++;
  169.         if(fwrite(obuf, uimg.pxl_out, fsizeNew, out_fp) != fsizeNew)
  170.             syserr("error during write h");
  171.     }
  172.     break;
  173. default:
  174. case 'd':{
  175. int    big_row, big_cln;
  176. byte    *tbuf;
  177.     message("diagnal: angle = %.2f parallel %c axle\n", angle, xyz);
  178.     angle *= M_PI / 180.;
  179.     cosl = cos(angle);
  180.     ssin = sin(angle);
  181.     sinl = fabs(ssin);
  182.     tanl = sinl / cosl;
  183.     if (bgn_frm > frms)    bgn_frm = frms;
  184.     if (bgn_row > rows)    bgn_row = rows;
  185.     if (bgn_cln > cols)    bgn_cln = cols;
  186.     switch(xyz){
  187.     case 'x':    /* max fsize = sqrt(r^2 + f^2) * c    */
  188.         big_row = sqrt((double)(rows*rows + frms*frms));
  189.         big_cln = cols;
  190.         break;
  191.     case 'y':
  192.         big_row = rows;
  193.         big_cln = sqrt((double)(frms*frms + cols*cols));
  194.         break;
  195.     case 'z':
  196.     default:
  197.         big_row = frms;
  198.         big_cln = sqrt((double)(rows*rows + cols*cols));
  199.     }
  200.     obuf = nzalloc(big_row, big_cln, "d_obuf");
  201.     if (frames>1){
  202.         tbuf = (byte*)zalloc(big_row, big_cln, "tbuf");
  203.         message("FramingR %d, FramingC %d\n", big_row, big_cln);
  204.     }
  205.     for (f=0; f<frames; f++){
  206.         switch(xyz) {
  207.         case 'x':
  208.         /*###############################################
  209.         #    column is same as old one        #
  210.         &    Z axle angle = 0. Parallel to X axle    #
  211.         #    Z parameter (frame #) is main parameter    #
  212.         ###############################################*/
  213.         width = cols;    l_r = fsize;
  214.         if (bf){
  215.             if (angle > 0){
  216.             rowsNew = (frms - bgn_frm + 1) / cosl;
  217.             bgn_row = 1;
  218.             }
  219.             else{
  220.             rowsNew = bgn_frm / cosl;
  221.             bgn_row = bgn_frm * tanl;
  222.             if (bgn_row > rows)
  223.             {
  224.                 bgn_frm = (bgn_row - rows + 1)/tanl;
  225.                 bgn_row = rows;
  226.             }
  227.             else    bgn_frm = 1;
  228.             }
  229.             if (rowsNew * sinl > rows)    rowsNew = rows / sinl;
  230.         }
  231.         else{    /* bgn_row mode    */
  232.             rowsNew = (rows - bgn_row + 1) / sinl;
  233.             if (rowsNew*cosl > frms)
  234.                 rowsNew = frms/cosl;
  235.             if (angle>0)    bgn_frm = 1;
  236.             else{
  237.                 bgn_frm = frms - (int)(rowsNew * cosl);
  238.                 bgn_row += rowsNew * sinl;
  239.             }
  240.         }
  241.         message("rn=%d, br=%d, bf=%d\n", rowsNew, bgn_row, bgn_frm);
  242.         clnsNew = cols;
  243.         obp = (byte*)obuf + clnsNew;
  244.         inbp = (byte*)inbuf+(bgn_frm-1)*fsize;
  245.         ibp = ibuf+fsize;
  246.         /* 1st coln is integer. There is no need of interpolating    */
  247.         memcpy(obuf, inbp+(bgn_row-1)*cols, clnsNew);
  248.         for (i=df=0; i<fsize; i++)    *ibp++ = *inbp++;
  249.         for(r = 1; r < rowsNew; r++)    /* rowsNew = frames    */
  250.         {
  251.             cfmod = r * cosl - df;
  252.             if (cfmod>=0)
  253.             {    ibp = ibuf+fsize;    df++;
  254.             memcpy(ibuf, ibp, fsize*sizeof(*ibp));
  255.             for(i = 0; i < fsize; i++)
  256.                 *ibp++ = (int)(*inbp++ & 0xFF);
  257.             }
  258.             else    cfmod += 1;
  259.             i = rfmod = r * ssin;
  260.             ibp = ibuf + (bgn_row-1 + i) * cols;
  261.             rfmod -= i;
  262.             for(c = 0; c < clnsNew; c++)
  263.             *obp++ = interpolation(ibp++, width, l_r, fabs(rfmod), cfmod);
  264.         }
  265.     break;
  266.         /*###############################################
  267.         #    new rows is old row on X-Y surface    #
  268.         &    Z axle angle = 0. For consistency.    #
  269.         #    new column is old row. Parallel Y axle    #
  270.         ###############################################*/
  271.         case 'y':
  272.         width = 1;    l_r = fsize;
  273.         if (bf){
  274.             if (angle > 0){
  275.             rowsNew = (frms - bgn_frm + 1) / cosl;
  276.             bgn_cln = 1;
  277.             }
  278.             else{
  279.             rowsNew = bgn_frm / cosl;
  280.             bgn_cln = bgn_frm * tanl;
  281.             if (bgn_cln > cols)
  282.                {    bgn_frm = (bgn_cln - cols + 1)/tanl;
  283.                 bgn_cln = cols;
  284.                }
  285.             else    bgn_frm = 1;
  286.             }
  287.             if (rowsNew*sinl>cols)    rowsNew = cols/sinl;
  288.         }
  289.         else{    /* bgn_cln mode    */
  290.             if (!bgn_cln)    bgn_cln++;
  291.             rowsNew = (cols - bgn_cln + 1) / sinl;
  292.             if (rowsNew*cosl > frms)    rowsNew = frms/cosl;
  293.             if (angle>0)    bgn_frm = 1;
  294.             else{
  295.             bgn_frm = frms - (int)(rowsNew * cosl);
  296.             bgn_cln += rowsNew * sinl;
  297.             }
  298.         }
  299.         message("nr=%d, bc=%d, bf=%d\n", rowsNew, bgn_cln, bgn_frm);
  300.         clnsNew = rows;
  301.         obp = (byte*)obuf;
  302.         inbp = (byte *)inbuf+(bgn_frm-1)*fsize;
  303.         ibp = ibuf+fsize;
  304.         for (i=0; i<fsize; i++)    *ibp++ = *inbp++;
  305.  
  306.         for(r=df=0; r < rowsNew; r++)
  307.         {
  308.             cfmod = r * cosl - df;
  309.             if (cfmod>=0){
  310.             ibp = ibuf+fsize;    df++;
  311.             memcpy(ibuf, ibp, fsize*sizeof(*ibp));
  312.             for(i = 0; i < fsize; i++)
  313.                 *ibp++ = (int)(*inbp++ & 0xFF);
  314.             }
  315.             else    cfmod += 1;
  316.             i = rfmod = r * ssin;
  317.             ibp = ibuf + bgn_cln-1 + i;
  318.             rfmod -= i;
  319.             for(c=0; c < clnsNew; c++, ibp+=cols)
  320.             *obp++ = interpolation(ibp, width, l_r, fabs(rfmod), cfmod);
  321.         }
  322.     break;
  323.         /*###############################################
  324.         #    new column is on the X-Y surface    #
  325.         &    X axle angle degree = 0.        #
  326.         #    new row is old frame. Parallel Z axle    #
  327.         ###############################################*/
  328.         default:
  329.         case 'z':    width = cols;
  330.         if (bc){
  331.             if (angle > 0){
  332.             clnsNew = (cols - bgn_cln + 1) / cosl;
  333.             bgn_row = 1;
  334.             }
  335.             else{
  336.             clnsNew = bgn_cln / cosl;
  337.             bgn_row = clnsNew * sinl;
  338.             if (bgn_row > rows)
  339.             {    bgn_cln = (bgn_row-rows+1) / tanl;
  340.                 bgn_row = rows;
  341.             }
  342.             else    bgn_cln = 1;
  343.             }
  344.             if (clnsNew * sinl > rows)    clnsNew = rows / sinl;
  345.         }
  346.         else{
  347.             clnsNew = (rows - bgn_row + 1) / sinl;
  348.             if (clnsNew * cosl > cols)    clnsNew = cols / cosl;
  349.             if (angle>0)    bgn_cln = 1;
  350.             else{
  351.             bgn_cln = cols - clnsNew * cosl;
  352.             bgn_row += clnsNew *sinl;
  353.             }
  354.         }
  355.         message("nc=%d, bc=%d, br=%d\n", clnsNew, bgn_cln, bgn_row);
  356.         rowsNew = frms;
  357.         obp = (byte*)obuf;
  358.         inbp = (byte*)inbuf;    /* convert to integer    */
  359.         for(r=0; r < rowsNew; r++)    /* rowsNew = frames    */
  360.         {
  361.             ibp = ibuf;
  362.             for(i=0; i < fsize; i++)
  363.             *ibp++ = (int)(*inbp++ & 0xFF);
  364.             ibp = ibuf + (bgn_row-1)*cols + (bgn_cln-1);
  365.             *obp++ = *ibp;
  366.             for(c=1; c < clnsNew; c++){    /* colsNew = ?/cosl */
  367.             register int    j = rfmod = c * ssin;
  368.             i = cfmod = c * cosl;
  369.             *obp++ = interpolation(ibp+j*cols+i, width, l_r,
  370.                 rfmod-j, cfmod-i);
  371.             }
  372.         }
  373.     }/* end of switch(xyz)    */
  374.  
  375.     if (!f){    /* first frame */
  376.     register int    r=rows, c=cols, f=frms;
  377.         if (frames > 1){
  378.             cols=big_cln;    rows=big_row;
  379.         }
  380.         else {    cols=clnsNew;    rows=rowsNew;    }
  381.         frms=frames;
  382.         fsizeNew = rows*cols;
  383.         (*uimg.header_handle)(HEADER_WRITE, &uimg, argc, argv, True);
  384.         frms=f;    rows=r;    cols=c;
  385.     }
  386.     if (bf)    bgn_frm += Sign(angle);
  387.     else if(bc)    bgn_cln += Sign(angle);
  388.     else    bgn_row += Sign(angle);
  389.     {
  390.     register byte*    opp=(byte*)obuf;
  391.     if (frames>1){
  392.         opp = tbuf;
  393.         framing(obuf, opp, rowsNew, clnsNew, big_row, big_cln);
  394.         if (bc && (!bgn_cln || bgn_cln > cols) ||
  395.             !bgn_row || bgn_row > rows ||
  396.             bf && (!bgn_frm || bgn_frm > frms))
  397.             f = frames;    /*    STOP -- out of boundary    */
  398.     }
  399.     i = fwrite(opp, uimg.pxl_out, fsizeNew, out_fp);
  400.     if(i!=fsizeNew)    syserr("error during write (%c)%d", xyz, i);
  401.     }
  402.     }
  403. }
  404. }
  405. exit(0);
  406. }
  407.  
  408.  
  409. /*=======================================================================
  410. *    cur_p is relative to current position.        |        *
  411. *    width:                        | 10 (-a) 11    *
  412. *        for X-Y (Z axle) is column - offset    |        *
  413. *        for Y-Z (X axle) is row - offset    | 00    01    *
  414. *        for Z-X (Y axle) is frames - ofset    |        *
  415. *    riF, ciF is exactly current position.        | 10    11    *
  416. *======================================================================*/
  417. interpolation(cur_p, span, l_r, uh, rw)
  418. register int    *cur_p;
  419. float    uh, rw;    /* the fraction percentage    */
  420. {
  421. register float    fh = 1.0 - uh,    /* the % of regular position    */
  422.         fw = 1.0 - rw,
  423.     p00, p01, p10, p11;
  424.  
  425. p00 = *cur_p * fh;
  426. p10 = *(cur_p + span) * uh;
  427. p01 = *(cur_p + l_r) * fh;
  428. p11 = *(cur_p + span + l_r) * uh;
  429. return    (fw*(p00 + p10) + rw*(p01 + p11));
  430. }
  431.  
  432. framing(ibp, obp, ir, ic, or, oc)
  433. byte    *ibp, *obp;
  434. {
  435. int    mtop = (or - ir) >> 1,
  436.     mleft= (oc - ic) >> 1;
  437. register int    r, c;
  438.  
  439. for (r=0; r<mtop; r++)
  440.     for (c=0; c<oc; c++)
  441.     *obp++ = 0;
  442. for (r=0; r<ir; r++){
  443.     for (c=0; c<mleft; c++)
  444.         *obp++ = 0;
  445.     for (c=0; c<ic; c++)
  446.         *obp++ = *ibp++;
  447.     for (c=0; c<oc-ic-mleft; c++)
  448.         *obp++ = 0;
  449. }
  450. for (r=0; r<or-ir-mtop; r++)
  451.     for (c=0; c<oc; c++)
  452.     *obp++ = 0;
  453. }
  454.